Home Range Analysis of Eastern Box Turtles

data <- read_csv("Eastern Box Turtle Tracking.csv")

data <- data[!is.na(data$'location-long'), ]
colnames(data)[colnames(data)=="location-long"] <- "long"
colnames(data)[colnames(data)=="location-lat"] <- "lat"
colnames(data)[colnames(data)=="individual-local-identifier"] <- "individual"
colnames(data)[colnames(data)=="utm-easting"] <- "easting"
colnames(data)[colnames(data)=="utm-northing"] <- "northing"

Eastern Box Turtle Movement in North Carolina

Data Collected by N.C. Wildlife Commission

Data accessed through www.movebank.org

Plot of Points

ggplot(data=data, aes(x = long, y=lat, color=individual))+ 
  geom_point()+
    labs(x="Longitude", y="Latitude") +
    guides(color=guide_legend("Identifier"))

lapply(split(data, data$individual), 
       function(x)write.csv(x, file = paste(x$individual[1],".csv"), row.names = FALSE))

files <- c("Billy Bob .csv", "Bon Jovi (BJ) .csv", "SAMMY .csv")

utm_points <- cbind(data$easting, data$northing)
utm_locations <- SpatialPoints(utm_points, 
                 proj4string=CRS("+proj=utm +zone=17 +datum=WGS84"))
proj_lat.lon <- as.data.frame(spTransform(
                utm_locations, CRS("+proj=longlat +datum=WGS84")))
colnames(proj_lat.lon) <- c("x","y")
raster <- openmap(c(max(proj_lat.lon$y)+0.001, min(proj_lat.lon$x)-0.001), 
                  c(min(proj_lat.lon$y)-0.001, max(proj_lat.lon$x)+0.001), 
                  type = "bing")

raster_utm <- openproj(raster, 
              projection = "+proj=utm +zone=17 +ellps=WGS84 +units=m +no_defs")

autoplot(raster_utm, expand = TRUE) + theme_bw() +
  theme(legend.position="bottom") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_point(data=data, aes(easting,northing,
             color=individual), size = 3, alpha = 0.8) +
  theme(axis.title = element_text(face="bold")) + labs(x="Easting",
        y="Northing") + guides(color=guide_legend("Identifier"))

Minimum Convex Polygon

mcp_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$easting)
  y <- as.data.frame(data$northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=17 +ellps=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  mcp.out <- mcp(xy, percent=100, unout="ha")
  mcp.points <- cbind((data.frame(xy)),data$individual)
  colnames(mcp.points) <- c("x","y", "identifier")
  mcp.poly <- fortify(mcp.out, region = "id")
  units <- grid.text(paste(round(mcp.out@data$area,2),"ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  mcp.plot <- autoplot(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=mcp.out, aes(x=mcp.poly$long, y=mcp.poly$lat), alpha=0.5, color = "red", fill = "red") +
    geom_point(data=mcp.points, aes(x=x, y=y)) + 
    labs(x="Easting (m)", y="Northing (m)", title=mcp.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  mcp.plot}


pblapply(files, mcp_raster)

  |                                                                                        | 0 % ~calculating  
  |==============================                                                          | 33% ~03s          
  |===========================================================                             | 67% ~01s          
  |========================================================================================| 100% elapsed=03s  
[[1]]

[[2]]

[[3]]

Kernel-Density Estimation

kde_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$easting)
  y <- as.data.frame(data$northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=17 +ellps=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=100)
  ver <- getverticeshr(kde, 95)
  ver1 <- getverticeshr(kde, 75)
  ver2 <- getverticeshr(kde, 50)
  kde.points <- cbind((data.frame(data.proj@coords)),data$individual)
  colnames(kde.points) <- c("x","y","identifier")
  kde.poly <- fortify(ver, region = "id")
  kde.poly1<- fortify(ver1, region = "id")
  kde.poly2<- fortify(ver2, region ="id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  kde.plot <- autoplot(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=kde.poly, aes(x=kde.poly$long, y=kde.poly$lat), alpha=0.5, color = "red", fill = "red") +
    geom_polygon(data=kde.poly1, aes(x=kde.poly1$long, y=kde.poly1$lat), alpha=0.5, color = "green", fill = "green") +
    geom_polygon(data=kde.poly2, aes(x=kde.poly2$long, y=kde.poly2$lat), alpha=0.5, color = "blue", fill = "blue") +
    geom_point(data=kde.points, aes(x=x, y=y)) +
    labs(x="Easting (m)", y="Northing (m)", title=kde.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  kde.plot
}

pblapply(files, kde_raster)

  |                                                                                        | 0 % ~calculating  
  |==============================                                                          | 33% ~01s          
  |===========================================================                             | 67% ~01s          
  |========================================================================================| 100% elapsed=02s  
[[1]]

[[2]]

[[3]]